% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication "Deconstructing the Yield Curve"
% Crump and Gospodinov (2024)
% Date: 26-JUL-2024
% Code for probability of recession application (Figure 6)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
addpath('functions/');
mex 'functions/buildVAR.c'
mex 'functions/buildForwards.c'

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OPTIONS
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
draft = 5;
pba = [2.5,1,1];
cvg = 0.68; 
dataSet = 'gsw';
freq = 'q';
N = 40;
B = 999;
h = 4;
doLegend = false;
startDate = 197103;
endDate = 202401;

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specifications
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
spec = 'standard';   % options are 'standard' or 'wright' or 'alt'
ytrans = @(x) x <= 0; ht = 'Single';
%ytrans = @(x) movmean((x <=0),[h-1,0])>0; ht = 'Cumulative';

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fixed parameters
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
l = 0;
p.M                    = 'ROT';
p.CG_B_BC              = 1000;

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOAD GDP DATA
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
rawGDP = readtable('../input/routput_first_second_third.xlsx','Sheet','DATA','Rang', 'A5');
gdpFinal = rawGDP.Third;
gdpFinal(isnan(gdpFinal)) = rawGDP.Most_Recent(isnan(gdpFinal));
%
yearGDP = [1965; 1965; kron((1966:round(endDate/100))',ones(4,1))];
qtrGDP = [3; 4; kron(ones(round(endDate/100)-1966+1,1),[1,2,3,4]')];
datesGDP = yearGDP*100+qtrGDP;
%
yearEndGDP = str2double(rawGDP.Date{end}(1:4));
qtrEndGDP = str2double(rawGDP.Date{end}(7));
%
datesGDP = datesGDP(datesGDP<=100*yearEndGDP+qtrEndGDP);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOAD BOND DATA
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
dataAll = loadBondData(N);
dataTag = ['dataAll.', dataSet, '.', freq];
for i = {'prc', 'yld', 'ret', 'xret', 'fwd', 'dret'}, eval([i{1}, 'Raw = ', dataTag, '.', i{1}, ';']);  end
data.mats =  eval([dataTag, '.mats(1:N)']);
datesBonds = eval([dataTag, '.dates']);
data.period = eval([dataTag, '.period']);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CLEAN DATA
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
dates = intersect(datesBonds, datesGDP);
if strcmp(dataSet,'gsw') && strcmp(freq,'q'), dates = dates(dates>=197103); end
idxStart = find(dates==startDate);
idxEnd = find(dates==endDate);
dates = dates(idxStart:idxEnd);
idxBonds = ismember(datesBonds, dates);
idxGDP = ismember(datesGDP, dates);

data.prc = prcRaw(idxBonds,1:N);
data.yld = yldRaw(idxBonds,1:N);
data.fwd = fwdRaw(idxBonds,1:N);
data.dret = dretRaw(idxBonds,1:N);
data.ret = retRaw(idxBonds,1:N);
data.xret = xretRaw(idxBonds,1:N);

data.pred = gdpFinal(idxGDP,:);
data.T = size(data.prc,1);
data.N = N;
data.h = h;

if strcmp(p.M,'ROT'), p.M=ceil((data.T*N)^(2/5)); end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATE PREDICTORS
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if strcmp(spec, 'standard')
    Xraw = data.yld(:,end)-data.yld(:,1);
elseif strcmp(spec, 'wright')
    Xraw = [data.yld(:,end)-data.yld(:,1), data.yld(:,1)];
elseif strcmp(spec, 'alt')    
    Xraw = [data.yld(:,end)-data.yld(:,1), data.yld(:,1)-movmean(data.yld(:,1),[11,0])];
end

yraw = ytrans(data.pred);

Xsample = Xraw(1:end-h,:);
XsampleFcst = Xraw;
ysample = yraw(h+1:end);

betaFullSample = glmfit(Xsample,ysample,'binomial','link','probit');
yfitFullSample = normcdf(betaFullSample(1) + XsampleFcst*betaFullSample(2:end));
betaBSdiffFullSample = NaN(numel(betaFullSample),B);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BOOTSTRAP DATA
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
bsYld = NaN([size(data.yld),B]);
bsGDP = NaN([numel(data.pred),B]);

for b = 1:B
    
    varJoint = var1BC([data.fwd(:,end) data.pred],p.CG_B_BC,0);
    Zfwd = varJoint.VBC(2:end,1);
    ZGDP = varJoint.VBC(2:end,2);

    Zinit = [Zfwd, data.dret(2:end,2:end), ZGDP];
    Z = [Zinit; Zinit(1:p.M-1,:)];
    idxRand = randi(data.T-p.M,[ceil(data.T/p.M),1]);
    idxBootstrap = kron(idxRand,ones(p.M,1)) + kron(ones(ceil(data.T/p.M),1),(0:p.M-1)');
    idxBootstrap = idxBootstrap(1:data.T-1);
    bsZ = Z(idxBootstrap,:);
    bsDret = [NaN(1,size(data.fwd,2)); NaN(size(data.fwd,1)-1,1) bsZ(:,2:N)];

    bsV = bsZ(:,[1,end])-ones(data.T-1,1)*mean(bsZ(:,[1,end]));
    
    bsJoint = buildVAR(varJoint.muBC,varJoint.PhiBC,bsV',[data.fwd(1,end) data.pred(1,:)]')';
    bsJoint = [[data.fwd(1,end) data.pred(1,:)]; bsJoint];

    bsFwdN = bsJoint(:,1);
    bsGDP(:,b) = bsJoint(:,2);

    bsFwd = buildForwards(bsDret, data.fwd(1,:), bsFwdN);
    bsYld(:,:,b) = cumsum(bsFwd,2)./cumsum(ones(size(bsFwd)),2);

    if strcmp(spec, 'standard')
        XrawBS = bsYld(:,end,b)-bsYld(:,1,b);
    elseif strcmp(spec, 'wright')
        XrawBS = [bsYld(:,end,b)-bsYld(:,1,b), bsYld(:,1,b)];
	elseif strcmp(spec, 'alt')
        XrawBS = [bsYld(:,end,b)-bsYld(:,1,b), bsYld(:,1,b)-movmean(bsYld(:,1,b), [11,0])];
    end
    yrawBS = ytrans(bsGDP(:,b));
    
    betaBSdiffFullSample(:,b) = glmfit(XrawBS(1:end-h+1,:),yrawBS(h:end),'binomial','link','probit')-betaFullSample;
    
end

betaFullSampleBC = betaFullSample - mean(betaBSdiffFullSample,2);
betaFullSampleBCupper = betaFullSample - quantile(betaBSdiffFullSample,(1-cvg)/2,2);
betaFullSampleBClower = betaFullSample - quantile(betaBSdiffFullSample,1-(1-cvg)/2,2);

beta = NaN(size(Xraw,2)+1,data.T-h);
betaBSdiff = NaN(size(Xraw,2)+1,data.T-h,B);
betaBC = NaN(size(Xraw,2)+1,data.T-h);
betaBCupper = NaN(size(Xraw,2)+1,data.T-h);
betaBClower = NaN(size(Xraw,2)+1,data.T-h);

yfitBC = NaN(size(yfitFullSample));
yfitBCupper = NaN(size(yfitFullSample));
yfitBClower = NaN(size(yfitFullSample));

if (l>=1), error('Code only allows for leave one out'); end

for s = 1:data.T-h
    if mod(s,5)==0, disp([num2str(s), ' of ', num2str(data.T-h-2*l)]); end
    Xtmp = Xsample;
    ytmp = ysample;
    Xtmp(s,:) = [];
    ytmp(s) = [];
    
    if (s==1)
        beta(:,s) = glmfit(Xtmp,ytmp,'binomial','link','probit');                
    else
        beta(:,s) = glmfit(Xtmp,ytmp,'binomial','link','probit','B0',beta0);                
    end    
    beta0 = beta(:,s);

    for b = 1:B
        if strcmp(spec, 'standard')
            XrawBS = bsYld(:,end,b)-bsYld(:,1,b);
        elseif strcmp(spec, 'wright')
            XrawBS = [bsYld(:,end,b)-bsYld(:,1,b), bsYld(:,1,b)];
        elseif strcmp(spec, 'alt')
            XrawBS = [bsYld(:,end,b)-bsYld(:,1,b), bsYld(:,1,b)-movmean(bsYld(:,1,b), [11,0])];
        end
        yrawBS = ytrans(bsGDP(:,b));

        XsampleBS = XrawBS(1:end-h,:);
        ysampleBS = yrawBS(h+1:end);        
        XtmpBS = XsampleBS;
        ytmpBS = ysampleBS;
        XtmpBS(s,:) = [];
        ytmpBS(s) = [];
                
        betaBSdiff(:,s,b) = glmfit(XtmpBS,ytmpBS,'binomial','link','probit')-beta(:,s);
    end 
    
    betaBC(:,s) = beta(:,s) - mean(betaBSdiff(:,s,:),3);
    betaBCupper(:,s) = beta(:,s) - quantile(betaBSdiff(:,s,:),(1-cvg)/2,3);
    betaBClower(:,s) = beta(:,s) - quantile(betaBSdiff(:,s,:),1-(1-cvg)/2,3);

    yfitBC(s) = normcdf(betaBC(1,s-l) + XsampleFcst(s,:)*betaBC(2:end,s-l));

    % Note: yfit and Xsample(s,:) are evaluated at s not s-l
    yfitBCupper(s) = normcdf(betaBCupper(1,s) + XsampleFcst(s,:)*betaBCupper(2:end,s));
    yfitBClower(s) = normcdf(betaBClower(1,s) + XsampleFcst(s,:)*betaBClower(2:end,s));
    
end

for s = data.T-h+1:data.T
    yfitBC(s) = normcdf(betaBC(1,data.T-h) + XsampleFcst(s,:)*betaBC(2:end,data.T-h));
    yfitBCupper(s) = normcdf(betaBCupper(1,data.T-h) + XsampleFcst(s,:)*betaBCupper(2:end,data.T-h));
    yfitBClower(s) = normcdf(betaBClower(1,data.T-h) + XsampleFcst(s,:)*betaBClower(2:end,data.T-h));
end

t1 = datetime('30-Jun-1971');
t2 = dateshift(t1,'end','quarter', (1:numel(dates)+h)');
datesPlot = t2(h+1:end);

save(['../output/probRec_Spec_', spec, '_Ht_', ht]);

figName = figure;
hold on;
h1 = plot(datesPlot, yfitFullSample, 'r', 'LineWidth', 1.5);
h2 = plot(datesPlot, yfitBC, 'k', 'LineWidth', 1.5);
pbaspect(pba);
set(gca, 'FontSize', 14);
recessionplot;
if doLegend, legend([h1,h2], {'MLE', 'Bias Corrected'}, 'FontSize', 14); end
hold off;
xlim([datetime('1-Mar-1971'), datetime('1-Dec-2025')]);
print(figName, '-dpdf', ['../output/probRec_', spec, '_transform', ht '_h', num2str(h), '_supp_draft', num2str(draft), '.pdf']);
close all;

figName = figure;
hold on
h1 = plot(datesPlot, yfitBCupper, '--k', 'LineWidth', 1.25);
h2 = plot(datesPlot, yfitBClower, '--k', 'LineWidth', 1.25);
tmpDates = [datesPlot', fliplr(datesPlot')];
tmpFit = [yfitBClower', fliplr(yfitBCupper')];
fill(tmpDates, tmpFit, [1,.8,.8]);
h1 = plot(datesPlot, yfitBC, 'k', 'LineWidth', 1.75);
pbaspect(pba);
set(gca, 'FontSize', 14);
recessionplot;
hold off;
xlim([datetime('1-Mar-1971'), datetime('1-Dec-2025')]);
print(figName, '-dpdf', ['../output/probRec_', spec, '_transform', ht '_h', num2str(h), '_draft', num2str(draft), '.pdf']);
close all;
